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AN ALGORITHM AND COMPUTER PROGRAM TO LOCATE REAL ZEROS 


OF REAL POLYNOMIALS 


David R . Hedgley , Jr . 
Flight Research Center 


INTRODUCTION 


Finding the real zeros of polynomials is a classical problem in almost every 
technical discipline (ref. 1) . It has assumed major importance in the last two 
decades in the treatment of the masses of data that have accompanied the growth of 
technology . 

Many methods for locating the zeros of real polynomials are being used. How- 
ever , all of these methods , which locate the zeros of a real polynomial with no prior 
information regarding the location of the zeros , find both the real and the complex 
zeros . Moreover , many of the methods have inherent weaknesses . For example , 
the polynomial x 20 - 1 causes the Newton-Raphson approach to diverge near 1 or -1. 
Bairstow's method requires close approximations to a zero; otherwise the results may 
be erroneous. Laguerre's method is satisfactory if the polynomial has all real zeros; 
however , if it has complex zeros , little can be said about its behavior . Reference 2 
discusses these methods in more detail. Finally, the Jenkins-Traub algorithm, which 
is considered the most advanced method , has difficulty with zeros which form a clus- 
ter . 


In addition to these anomalies , the methods are inefficient when only real zeros 
are desired. Furthermore, because of possible computational inaccuracies, which 
are to a large degree a function of word size limitations of computers , real zeros 
can be mistaken for complex zeros when the imaginary part of the zero is small in 
absolute value. 

The intent of this paper is to present an algorithm which (1) presupposes no 
knowledge of the location or number of real zeros and (2) compares favorably with 
the standard methods when a polynomial has all real zeros but (3) demonstrates a 
pronounced superiority in efficiency when the polynomial has complex zeros . 

A computer program to implement the algorithm is presented , and results from 
the Laguerre, Newton-Raphson, and Jenkins-Traub methods are compared with 
results obtained from the proposed method . 



BACKGROUND 


Three significant criteria for evaluating a technique which locates real zeros 
of a real polynomial are: its inherent rate of convergence, the computational time 
required for each iteration , and the probability of convergence . The two-point 
secant method for locating real zeros is given by the equation 



in which x ^ is an arbitrary variable on the real axis, x^ is the next iterate deter- 
mined by x j and x^, and p^x^ is value of a real polynomial, p(x) , at x„ This 

method has an excellent rate of convergence, (1 + V5)/2 (ref. 2). The computational 
time per iteration is small because only P(^ +1 ) need be computed after the first 

iteration . Furthermore , this formulation of the secant method enhances numerical 
stability , since only a few significant digits are required as convergence is neared 
(ref. 2) . Therefore, because this method satisfies the convergence and computa- 
tional time criteria, it was selected for further development of the proposed algo- 
rithm . 

The following sketch is a geometrical representation of the iterative process of 
the two-point secant method: 


y 



In the equation for the secant method, initially some estimates for a\ and x fl are 
made. Geometrically, x^ is the x-intercept of a straight line determined by 
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(V( oc f ^J and |“ < 3c-_ 1 .p^jJJ , which is a secant through the curve p (x) . Having 
determined ^ +1 » x^ and x f+ ^ are used to generate x. +2> and so on. 


Although this process is efficient when it converges , it may not a: 


to a zero. If the curve is peculiar, it may diverge, that is, 



ways converge 
> 1 . 


If the zero, r, is simple, the secant method will converge to r for and x ^ 

sufficiently "close" to r, that is, < t, where t is an arbitrary tolerance 

limit (ref. 2) . However, if r is of multiplicity greater than 1, it may not converge 
for any choice of step size, A such that - r| < A^. and \x._^ - r | < A^. In 

fact, it may happen that for a\ ^ x i~i > P( x i 

leads to divergence. 


^ This condition obviously 


On the basis of this analysis of the strengths and weaknesses of the secant 
method , the proposed algorithm was developed , as discussed in the next section . 


ALGORITHM DEVELOPMENT 


Increasing the probability of convergence of the secant method requires a 
theorem on the bounds of real zeros of a real polynomial , p (oc) . The following 
theorem was adapted from reference 3: 


Theorem 1. If p(oc) = C n x n + C n _^x n ^ . . . + Cq, ^C n > 0^ is a real polyno- 
mial (where C^,i = 0 ,n are the coefficients of p(ac)) and if the first negative coeffi- 
cients are preceded by k coefficients which are positive or zero, and if g denotes 
the greatest number in absolute value of the negative coefficients , then each positive 

real zero is less than the quantity 1 + kj g/C 

This theorem makes it possible to find an interval, I, which contains all the 
real zeros of p Gx) , for , by the theorem , the lower bound for the real zeros of p (jc) 
is the negative of the upper bound for p (-cc) . 


We now define A^. in the following way: 


A . = 4 j = 1 , 2 , 3 , 
J 2 ] 


n 


where n = greatest integer less than or equal to log 2 I/e in which s is an arbitrarily 
small constant less than 1 , and I = max( I UB I , I LB I ) in which max is the larger of 
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the upper bound , UB , or the lower bound , LB , value of the real zeros of p (a?) . The 
± indicates that the iteration proceeds within the bounded interval in both the posi- 
tive and the negative direction and in an alternating manner . This scheme usually 
permits the extraction of real zeros in increasing magnitude , thus preserving 
accuracy (ref. 2) . 

Further, we define x ^ and x^ , the initial choices, with step size as follows: 


X i = LA j 


x._ t = (L - 1)A. 


L = 1,2, 3, 


, k; k = 2^ 


where j is held fixed . 

The process is completed if for A .,j = 1 and L = 1, x. and x^ lead to converg- 
ence. If not, L = 2 is selected, keeping A fixed, and the process is repeated. If, 
however , no choice of L for j = 1 leads to convergence , step size A^ is chosen where 

Ag = A^/2 by the preceding definition, and the previous steps are repeated. That 
is, every L for each A.,j = 1,2 n has the potential for convergence where 

n n-l n-Z l 


This iterative scheme for choosing initial values x ^ and x^_^ increases the 

probability that the secant method will converge, provided p(x ) has a real zero, 
for we know that the secant method will converge to r , a simple zero , for x . , and 

x._^ close to r. In fact, it may converge when r is not simple and x^ and x^ are 
not close to r . Clearly , if x ^ and x^ are assured of being close to r by decreasing 
A^., the probability of convergence is greater because of the large number of avail- 
able choices as well as the fact that at least one pair will be close to a real zero. 

Finally, if this scheme is not successful and if the existence of at least one real 
zero is assumed, it is highly probable that r is not a simple zero. Then consider 
the following reasoning . Let e be such that when Ay < e the secant method is no 

longer considered fruitful. If for Ay < e, p(x^ does not converge, find b such 

that j p(b)j = min|p^x^| , where min is the minimum value and x ^ is any value 

chosen or computed for which Ay < e and divergence occurred for every Ay > s . 

From the step size and because of the continuity of p (x) , we shall assume that a 
real zero , r , exists such that I b - r I < Ay . Consider the closed interval 

| 'b - A.,b + AyJ. Using b as a center and subjecting this interval to the bisection 

method whose direction of seek is governed by the absolute value of the end points 
at every subsequent subinterval, it is again probable that p(x.\ will converge to 
P(r). V l/ 
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Although the combination of the preceding schemes is efficient in locating a real 
zero, its efficiency and validity in determining if all the real zeros have been 
located or if any real zeros exist can be poor , particularly in higher degree poly- 
nomials with complex zeros. Sturm’s theorem (ref. 2) removes this difficulty. The 
theorem is stated as follows: 

Theorem 2. Let jp^(l ^ i < n + 1) be a sequence of polynomials related to 
p(x) of degree n in the following way: 

P 1 (oc) = p(x) P 2 C*) = P'(x ) 

P t _ iGxO = q i _ 1 (x)p. (x) - P i+1 (x) i=2,3,...,m-l<n 

Pm-l (5c) = P m -l to> Pm W 

where q^_^(x) is the quotient when p^Cxr) is divided by p f (x) , and p f+1 (x) is the 

negative of the remainder. If [e,f] is any interval on the real axis such that p(e) 
and p(f) f 0, then v(e) - v(f) is the number of distinct real zeros in [e,f] where 
v(e) and v(f) represent the number of variations of sign of jp^(x)J evaluated at e 
and f, respectively . 

Since at x = ±1, p(x) does not vanish, the implication is that the number of 
distinct real zeros of p(x) can be determined by using Theorem 2. Moreover, once 
a zero is located and the polynomial is deflated to give, for example, d(x) , a new 
I is determined using Theorem 1 and, hence, the number of distinct real zeros, if 
any, for d(x) can be determined, and so forth. 

Thus, for any real polynomial, p(x) , the status of completion with respect to 
all the real zeros including multiplicities can be ascertained efficiently and accu- 
rately with Theorems 1 and 2 by considering subsequent deflated polynomials and 
their corresponding interval of bounds in the same way . 


IMPLEMENTATION AND RESULTS 


The proposed algorithm was implemented by using an assembly of computer 
programs . Listings for the programs , together with brief descriptions and flow 
charts, are presented in the appendix. The programs were run on a Control Data 
Corporation Cyber 73-28 computer. The algorithm was applied to five polynomials: 
a fourth-order polynomial with a non-simple zero; a fourth-order polynomial with 
both simple and non-simple zeros; a thirteenth-order polynomial with all distinct 
real zeros which form a cluster; a fifteenth-order polynomial with complex zeros; and 
a twenty-fifth-order polynomial with complex zeros . The results are compared in 
tables 1 to 5 with results obtained on the CDC Cyber 73-28 computer for the Laguerre, 
Newton -Raph son , and Jenkins-Traub methods . A subroutine called ZPOLYR , 
obtained from International Mathematical and Statistical Laboratories , Inc . , was 
used to implement Laguerre's method. The Newton -Raphson method was imple- 
mented by using a subroutine called POLRT from the IBM Scientific Subroutine 
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Package . The Jenkins-Traub method was implemented by using the subroutine 
ZRPOLY , also taken from the International Mathematical and Statistical Laborato- 
ries , Inc . The zeros located by each of the methods are listed in the order found 


TABLE 1. -COMPARISON OF RESULTS OBTAINED FOR A FOURTH-ORDER POLYNOMIAL 
WITH A NON-SIMPLE ZERO 


Polynomial coefficients: 

. 10000E+01 
- , 40000E+01 
.60000E+01 
- . 40000E+01 
. 1 0000E+01 


Laguerre method 

Newton-Raphson method 

Jenkins-Traub method 

Proposed method 

Zeros 

Zeros 

Zeros 









Real zeros 

Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 


. 10000E+01 

. 80831E-07 

. 99998E+00 

0. 

. 10000E+01 

0. 

. 10000E+01 

. 10000E+01 

- . 80831E-07 

. 10000E+01 

0. 

. 10000E+01 

0. 

. 99996E+00 

. 10000E+01 

0. 

. 10000E + 01 

0. 

. 10000E+01 

0. 

. 10000E+01 

. 10000E+01 

0. 

. 10000E+01 

0. 

. 10000E+01 

0. 

.99998E+00 




Execution time, sec 



0.030 

0.195 

0. 

030 

0.064 


TABLE 2. -COMPARISON OF RESULTS OBTAINED FOR A FOURTH -ORDER POLYNOMIAL 

WITH ALL REAL ZEROS 


Polynomial coefficients: 

. 58500E+02 
- . 18000E+02 
- . 17500E+02 
. 40000E+01 
. 10000E+01 


Laguerre method 

Newton-Raphson method 

Jenkins-Traub method 

Proposed method 

Zeros 

Zeros 

Zeros 








Real zeros 

Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 


.21231E+01 

0. 

. 21231E+01 

0. 

. 21008E+01 

0. 

- . 21213E+01 

- . 61231E+01 

0. 

- . 61231E+01 

0. 

- . 21213E+01 

0. 

.21231E+01 

- . 21213E+01 

0. 

- . 21213E+01 

0. 

.21436E+01 

0. 

. 21213E+01 

.21213E+01 

0. 

. 21213E+01 

0. 

- . 61231E+01 

0. 

- . 61231E+01 



Execution time, sec 


0.042 

0.089 

0.031 

0.047 
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TABLE 3. -COMPARISON OF RESULTS FOR A THIRTEENTH-ORDER POLYNOMIAL 
WITH ALL DISTINCT REAL ZEROS WHICH FORM A CLUSTER 


Polynomial coefficients: 

. 30974E+03 
. 26695E+04 
. 10563E+05 
. 25410E+05 
. 41464E+05 
. 48468E+05 
. 41758E+05 
. 26849E+05 
. 12884E+05 
.45580E+04 
. 1 1554E+04 
. 19877E+03 
.20800E+02 
. 10000E+01 


Laguerre method 

Newton-Raphson method 

Jenkins-Traub method 

Proposed method 

Zeros 

Zeros 

Zeros 








■o * | m 








Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 


- . 20997E+0 1 

0. 

- . 10000E+01 

0. 

- . 1 1005E+01 

0. 

-.10G00E+01 

- . 22000E+01 

0. 

- . 11000E+01 

0. 

- . 99996E+00 

0. 

- . 11000E+01 

- . 20011E+01 

0. 

- . 12001E+01 

0. 

- . 13O91E+01 

0. 

- . 12000E+01 

- . 18970E+01 

0. 

- . 12995E+01 

0. 

- . 1 1976E+01 

0. 

- . 13000E+01 

- .18051E+01 

0. 

- . 14017E+01 

0. 

- . 15510E+01 

0. 

- . 14000E+01 

- . 16936E+01 

0. 

- . 14965E+01 

0. 

- . 15463E+01 

0. 

- . 15003E+01 

- . 16057E+01 

0. 

- . 16057E+01 

0. 

- . 13829E+01 

0. 

- . 15995E+01 

“ . 14965E+01 

0. 

- . 16936E+01 

0. 

- . 17802E+01 

0. 

- . 17005E+01 

- . 14017E+01 

0. 

- . 18051E+OI 

0. 

-.19065E+01 

0. 

- . 17997E+01 

- . 12995E+01 

0. 

- . 18970E+01 

0. 

- . 17277E+01 

0. 

-.19000E+01 

-.12001E+01 

0. 

- .20011E+01 

0. 

- . 21003E+01 

0. 

- , 20001E+01 

- . 11000E+01 

0. 

- . 20997E+01 

0. 

- . 19980E+01 

0. 

- . 21000E+01 

- . 10000E+01 

0. 

- .22000E+01 

0. 

- . 22000E+01 

0. 

- .22000E+01 

Execution time, sec 

0 

220 

0.950 

0 

285 

0.400 
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TABLE 4. -COMP ARISON OF RESULTS FOR A FIFTEENTH-ORDER POLYNOMIAL 
WITH COMPLEX ZEROS 


Polynomial coefficients: 

- . 10000E+01 
.20000E+01 
- . 30000E+01 
. 40000E+01 
- . 50000E+01 
. 60000E+01 
- . 70000E+01 
. 80000E+01 
-.90000E+01 
. 10000E+02 
- .11000E+02 
. 12000E+02 
- . 13000E+02 
. 14000E+02 
- . 15000E+02 
. 16000E+02 


Laguerre method 

Newton-Raphson method 

Jenkins-Traub method 

Proposed method 

Zeros 

Zeros 

Zeros 








T? - 








Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 


- . 56747E+00 

. 63594E+00 

.80860E+00 

0. 

. 80860E+00 

0. 

.80860E+00 

- . 56747E+00 

- . 63594E+00 

. 75103E+00 

- . 30205E+00 

- . 56747E+00 

. 63594E+00 


. 58601E+00 

- . 56241E+00 

. 75103E+00 

. 30205E+00 

- . 56747E+00 

- . 63594E+00 


. 58601E+00 

. 56241E+00 

- . 79267E+00 

- . 38555E+00 

- .28130E+00 

. 78671E+00 


.80860E+00 

0. 

- . 79267E+00 

.38555E+00 

- . 28130E+00 

- . 78671E+00 


- .28130E+00 

- . 78671E+00 

.33562E+00 

- . 74495E+00 

. 33562E+00 

. 74495E+00 


-.28130E+00 

- . 78671E+00 

. 33562E+00 

.74495E+00 

. 33562E+00 

- . 74495E+00 


- , 79267E+00 

- . 38555E+00 

. 33212E-01 

- . 82381E+00 

. 33212E-01 

. 82381E+00 


- . 79267E+00 

. 38555E+00 

. 33212E-01 

.82381E+00 

. 33212E-01 

- . 82381E+00 


. 33212E-01 

82381E+00 

- . 56747E+00 

- .63594E+00 

- . 79267E+00 

. 38555E+00 


.332I2E-01 

. 82381E+00 

- . 56747E+00 

.63594E+00 

- . 79267E+00 

- . 38555E+00 


. 75103E+00 

- .30205E+00 

. 58601E+00 

- .56241E+00 

. 58601E+00 

.56241E+00 


. 75103E+00 

.30205E+00 

. 58601E+00 

. 56241E+00 

. 58601E+00 

- . 56241E+00 


. 33562E+00 

. 74495E+00 

- .28130E+00 

- . 78671E+00 

. 75103E+00 

•30205E+00 


. 33562E+00 

74495E+00 

- . 28130E+00 

. 78671E+00 

. 75103E+00 

- . 30205E+00 


Execution time, sec 

19. 

910 

0.466 

0.208 

0.068 
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TABLE 5. -COMPARISON OF RESULTS FOR A TWENTY -FIFTH -ORDER POLYNOMIAL 
WITH COMPLEX ZEROS 


Polynomial coefficients: 


. 10000E+01 

. 14000E+02 

. 20000E+01 

. 15QQQE+02 

. 30000E+01 

. 16000E+02 

. 40000E+01 

. 17000E+02 

- . 44000E+02 

. 18000E+02 

. 60000E+01 

. 19000E+02 

.70000E+01 

. 20000E+02 

. 80000E+01 

. 21000E+02 

. 90000E+01 

. 22000E+02 

. 10000E+02 

. 23000E+02 

. 11000E+02 

. 24000E+02 

. 12000E+02 

.25000E+02 

. 13000E+02 

. 26000E+02 


Laguerre method 

Newton -Ra 

phson method 
ros 

Jenkins-Traub method 

Proposed method 

Zeros 

Ze 

Zeros 

Real zeros 

Real 

Imaginary 

Real 

Imaginary 

Real 

Imaginary 

No solution 

No solution 

- . 32799E+00 

0. 

-.57462E-01 

. 35331E+00 

- . 32799E+00 



- . 57462E-01 

- . 35331E+00 

- . 57462E-01 

- . 3533IE+00 

. 59100E+00 



- . 57462E-01 

. 35331E+00 

- . 32799E+00 

0. 

. 73831E+00 



. 59100E+00 

0. 

, 59100E+00 

0. 




. 73831E+00 

0. 

. 79772E+00 

. 61254E+00 




- . 10511E+01 

- . I5500E+00 

, 79772E+00 

- . 61254E+00 




- . 1051 1E+01 

. 15500E+00 

. 73831E+00 

0. 




. 33707E+00 

- . 98054E+00 

- , 54917E+00 

. 90186E+00 




. 33707E+00 

. 98054E+00 

- . 54917E+00 

- . 90186E+00 




. 91694E+00 

- . 33985E+00 

. 40625E-01 

. 10441E+01 




. 91694E+00 

. 33985E+00 

, 40625E-01 

- . 1Q441E+01 




- . 54917E+00 

- . 90186E+00 

. 59772E+00 

. 83228E+00 




- . 54917E+00 

. 90186E+00 

. 59772E+00 

- . 83228E+00 




- . 96101E+00 

- . 45127E+00 

. 33707E+00 

. 98054E+00 




- . 96101E+00 

. 45127E+00 

. 33707E+00 

- . 98054E+00 




- .26412E+00 

- . 10172E+G1 

-.78862E+00 

.70767E+00 




- . 26412E+00 

. 1Q172E+01 

- . 78862E+00 

- . 7O767E+O0 




. 59772E+00 

- . 83228E+00 

- . 26412E+00 

. 10172E+01 




. 59772E+00 

. 83228E+00 

- . 26412E+00 

- . 10172E+01 




. 79772E+00 

- . 61254E+00 

. 91694E+00 

. 33985E+00 




.79772E+00 

.61254E+00 

.91694E+00 

- . 33985E+00 




- . 78862E+00 

- . 70767E+00 

- . 96101E+00 

. 45127E+00 




- . 78862E+00 

. 70767E+00 

- . 96101E+00 

- . 45127E+00 




. 40625E-01 

- . 10441E+01 

- . 1051 1E+01 

. 15500E+00 




. 40625E-01 

. 10441E+01 

- . 10511E+01 

- . 15500E+00 



Execution time, sec 


1.559 J 0.438 0.238 


These results show that the proposed method compares favorably with the 
Laguerre, Newton-Raphson, and Jenkins-Traub methods when the polynomial has 
all real zeros , and is more efficient when the polynomial has complex zeros . More- 
over, as shown in table 1, Laguerre’s method identifies complex zeros when in 
fact the zeros are real. This discrepancy is possible with any method that locates 
real and complex zeros, thus demonstrating the advantage of a method which locates 
only real zeros . 
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EVALUATION OF ALGORITHM 


The proposed algorithm, by using Theorem 1 and Theorem 2 in combination, 
significantly reduces the difficulty of determining the number of real zeros of a 
polynomial and, hence, the status of completion. Additionally, the modified bisec- 
tion method , which is used when the secant method used iteratively does not lead 
to convergence for non-simple zeros , appreciably improves the probability of con- 
vergence . Since the predominant method is the secant method which has near 
quadratic convergence and small computational time per iteration , the proposed 
algorithm satisfies the previously stated criteria on rate and probability of converg- 
ence and computational time. 


Flight Research Center 

National Aeronautics and Space Administration 
Edwards, Calif., February 24, 1975 
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APPENDIX - COMPUTER PROGRAM 


PROGRAM DESCRIPTION 


The digital computer programs which implement the proposed algorithm were 
written in FORTRAN IV and occupy approximately 1000 decimal words excluding 
FORTRAN system routines. A user-written program calls the REALRT subroutine, 
which calls the remaining subroutines . 

SAMPLE j6 F USER-WRITTEN CALLING PROGRAM 

The following program , an example of a user-written program , constructs the 
polynomial x 4 - 8000or + 3x 2 + lOx 30 000 = p (x) and computes all the real 
zeros of this polynomial . 


5 


10 


15 


20 


25 


PROGRAM PET< INPUT, OUTPUT, TAPE3= OUTPUT, TAP E1=INPUT) 
(DIMENSION A (301 ,B(30) ,FOCT(30> 

DIMENSION 0(30) 

22 FORMAT ( 2F12, 4) 

A (1)=-30Q00 
A (21=10 

-A (3 ) =±- A^-3) S 


A (4)=-8000 
A 1 5 ) — 1 


N=4 

CALL REALRT (A, N, ROOT, NI) 

WRITE (3,7550) 

7550 FCRMATUH1) 

NN=N*1 
WRITE (3,24) 

DO 191 J=1,NN 
191 WRITE ( 3 ,76 ) A (J) 

WRITE (3,25) 

DO 190 J= 1 * N I 
WRITE (3,77) ROOT! J) 

190 CONTINUE 

24 FCRMAT( 10X»23H POLYNOMIAL COEFFICIENTS/) 

25 FORMAT (/ / 1 OX , 2 OH REAL ZEROES FOR PCX)/) 
78 FORMAT (3X ,E15. 5) 

77 FORMAT C3X,E15,5) 

STOP 

END 
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APPENDIX — Continued 


SUBROUTINES 


Subroutine REALRT (A, N, ROOT, NI) 

Purpose: To locate all the real zeros of a real polynomial of degree less than 30 
and greater than 1 

Flow chart: 



Subroutine arguments: 

A array of coefficients of p ( x ) 

N degree of p (oc) 

ROOT array containing the real zeros 

NI number of real zeros in array root 


Subroutine listing: 


5 


10 


15 


30 


SUBROUTINE RE ALRT ( A , N, ROOT , N I > 

C 

C SUBROUTINE REALRT LOCATES ALL REAL ZEROES OF A REAL 

C POLYNOMIAL WHOSE DEGREE IS LESS THAN 30 AND GREATER THAN 1. 

C 

C THE * A * ARRAY IS THE ARRAY CF COEFFICIENTS ARRANGED IN 

C ASCENOING CROER. 

C N IS THE DEGREE OF THE POLYNOMIAL 

C THE ROCT VARIABLE IS THE ARRAY WHICH WILL CONTAIN THE REAL 

C ROOTS. 

C NI INDICATES THE NUMBER OF REAL ZEROES FOUND. 

C 

DIMENSION A(1),RC0T(1) *9(30) 

NI= 0 

CALL BCUNE tA,N, BOUND, B > 

CALL CHECK ( A, N , BOUNO , I ER > 

IF( IER.EQ. 0) GO TO 120, 

N 1= IER 

CALL ROOEC A, B,N, BOUND, ROOT, NI) 

120 RETURN 
C 

END 
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APPENDIX - Continued 


Subroutine BOUNE (A, N, BOUND, B) 

Purpose: To determine an interval on the real axis that contains all the real zeros 
of p ( x ) 


Flow chart: 



Subroutine arguments: 


A 

N 

BOUND 

B 


array of coefficients for p (oc) 
degree of p Cxr) 

closed interval [-BOUND .BOUND] which contains all the real 
zeros of p (oc) , an arbitrary polynomial 
work array 
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APPENDIX - Continued 


Subroutine listing: 


5 


10 


15 


20 


25 


30 


35 


40 


45 


50 


55 


60 


65 


70 


SUBROUTINE 80UNE< AfN, BOUND, 9) 

c 

C THIS ROUTINE O ET ER HINES THE INTERVAL CONTAINING THE ROOTS, 

C IF ANY. {REFER TO CORRESPCNINC " TNX**, ) 

C 

DIMENSION E (1 ) , A 1 1 ) 

HHH=1Q.**6 - T * '' / ‘ ,c ' r 

K=N *1 

OO 15 J*i*K 

15 B(J)=A(J) 

C 

C FIND F { -X ) AND SUBSEQUENTLY COUNT THE NO. OF COEFFICIENTS 

C PRECEDING THE FIRST NEG COEFFICIENT. 

C 

C THEN 00 LIKEWISE FOR F<X>. 

C 

00 16 J=2,K 
T 1= (MOD ( J , 2 ) ) - • 5 

16 B(J>=B CJ I * <SIGNC1.,TL> ) 

C 

C IF LEADING COEFFICIENT IS NEG, THEN CONSIDER -F(X). 

C 

IF ( B (K > .GT *0)GO TO 190 
DO 22 J=1 , K- 
22 B(J)=-B(J) 

190 CONTINUE 
100 1=0 

00 2000 J=1 , K 

IF(9(K-J+1) .LT.O) GO TO 1200 

1 = 1*1 

20 0 G CONTINUE 
C 

C COMPUTE MAXIMUM NEG COEFFICIENT OF F(-X). 

C 

1200 R=0 

DO 1800 J=1,K 

IF(S( J).GE.O)GO TO 1800 

S = A BS ( B < J ) ) 

R=AMAX 1 C R, S) 

1800 CONTINUE 
W=l./I 

BT= 1.+ <R/B (K) 

C 

C COMPUTE MAXIMUM NEG COEFFICIENT CF F(X) 

C 

T = 0 

IF ( A (K ) . GT *0)GO TO 10 
T = 1 

DO 19 J= 1 » K 

19 A ( J) = ~A( J) 

10 1 = 0 

DO 20 J= 1, K 

IFtA (K-J +1 ) .LT • 0 l GO TO 12 
1 = 1+1 

20 CONTINUE 
12 R=Q 

DO 180 J=i,K 

IF ( A ( J) , GE • 0 ) GO TO 180 

S= ABS ( A ( J ) ) 

R = A MAXI (R, S) 

180 CONTINUE 
W=l./I 

BOUNO= 1. + ( { R/ A (K) ) **W> 

C 

C DETERMINE THE LARGER OF THE 8CINCS. THIS WILL BE THE INTERVAL 

C CONTAINING ROOTS. 

C 

BOUNO=AMAXl(BOUND,BT) 

IF(T.NE.l)GO TO 200 
00 11 J=l, K 
11 A ( J) = -A ( J ) 

200 RETURN 
END 
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APPENDIX — Continued 


BOUND 

IER 


closed interval [-BOUND .BOUND] which contains all the real zeros 
of p (x) , an arbitrary polynomial 
total number of distinct real zeros of p (x) 


Subroutine listing: 


5 


10 


15 


20 


25 


30 


35 


40 


45 


50 


55 


60 


65 


SUBROUTINE CHECK C A ,N . BOUNO • IER ) 

C 

C THIS ROUTINE C0NSTUCTS THE STURM SEQUENCE CF POLYNOMIALS 

C HAVING DETERMINED THIS SEQUENCE, IT THEN EVALUATES EACH 

C POLYNOMIAL AT BOUNO AN0 -BOUND ANC SUBSEQUENTLY COUNTS THE NUMBER OF 

C VARIATIONS OF SIGNS AT EACH EXTREMITY. THE DIFFERENCE 

C OF THESE SUMS REPRESENTS THE NUMBER OF REAL 0ISTINCT ZEROES 

C OF P(X). (I.E. V(A)-VIBI) 

DIMENSION AC1) ,B(30),C (30) ,0(30) 

IJ = 0 
K=N*1 
NN = K- 1 

XX=1./ABS (A (1) ) 

DO 22 JX=1,K 
22 0<JX)=A< JX)*XX 
FF=.1**10 
HH=10.**10 
DO 44 JX= 1 ,K 
GG= ABS (O C J X > ) 

HH= AMIN1 ( HH, GG ) 

44 CCNTINUE 
EPS = HH/4 

EPS=AMIN1(EPS,FF) 

CALL PDER(BfNN»A»K) 

1=0 
J = Q 

BO=BOUND 
0I = -8OINO 

CALL PVAL(S,BC,A,K) 

CALL PVAL(SO,BI,A,K) 

S=SIGN (l.,S) 

SC= SIGN (1 • ,S0 ) 

12 CALL PClVCC,NI,D*K,e,NN,SPS,IET> 

C 

C IF THERE IS NO REMAINDER* THEN COMPLETE ARGUMENT AND EXIT 

C 

IF (K.LE. 0 ) IJ=1 
CALL PVAL (T,80,8,NN) 

CALL PVAL (TC,EI,0,NN) 

T=S IGN ( 1 • , T) 

TO= SIG N ( 1 • , TO ) 

S=S*T 

SO=SO*TO 

IF (S. GT. 0. ) GO TO 20 
1=1*1 

20 IFtSO.GT. 0. )GO TO 21 
J=J*1 

21 CONTINUE 

IF< IJ.EQ.l )GO TO 10 

SC= TO 

S=T 

DC 17 JX= 1 , NN 

17 C(JX) = B(JX) 

DO 18 JX= 1 , K 

18 B(JX)=-D(JX) 

DO 19 JX= 1 , NN 

19 D(JXI=C(JX> 

L=K 

K=NN 

NN=L 

GO TO 12 
10 NG=J-I 
IER=NG 
RETURN 
END 
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APPENDIX - Continued 


Subroutine PDER (Y, IDIMY , X, IDIMX) 


Purpose: To find the derivative of a polynomial 


Subroutine arguments: 

Y array of coefficients , ordered from smallest to largest power for 

the derivative 
IDIMY dimension of Y 

X array of coefficients for the original polynomial 

IDIMX dimension of X 


Subroutine listing: 


SUBROUTINE POER ( Y f IOIM Y, X, IDIMX ) 
OIMENSION Xli),Y<l> 

IF(IDIMX-l) 3,3,1 

1 EXPT=0 

5 DO 2 1=1 , ICIMY 

EXPT=EXPT+1. 

2 Y <I) = X (1+1 ) *EXPT 
GO TO <♦ 

3 1 01 MY = 0 

10 <♦ RETURN 

ENO 


Subroutine PVAL (RES, ARG, X, IDIMX) 


Purpose: To evaluate a given polynomial at any given value using nested arithmetic 


Subroutine arguments: 


RES resultant value (i.e., P(ARG)) 

ARG given value of the independent variable 

X vector of coefficients , ordered from smallest to largest 

IDIMX degree + 1 

Subroutine listing: 


5 


10 


15 


SUBROUTINE PVfl L (RES , ARG , X, I DIMX ) 

C 

C THIS SUBROUTINE EVALUATES A GIVEN POLYNOMIAL AT 

C ANY VALUE USING NESTED ARITHMETIC 

C 
C 

OIMENSION X(i) 

RES = 0 
J= 1 0 IMX 

1 IF ( 3,3,2 

2 RES=RES*ARG+X< J) 

J = J-l 

GO TO 1 

3 RETURN 
END 
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APPENDIX - Continued 


Subroutine PDIV (P, IDIMP , X, IDIMX , Y, IDIM&, TOL, IER) 

Purpose: To divide two polynomials 

Subroutine arguments: 

P resultant vector of integral part 

IDIMP dimension of P 

X vector of coefficients for dividend polynomial, ordered from 

smallest to largest; replaced with remainder after division 
IDIMX dimension of X 

Y vector of coefficients of divisor polynomial ordered from smallest 

to largest 

IDIMY dimension of Y 

TOL tolerance value below which coefficients are eliminated 

IER error code; 1 denotes zero divisor , 0 denotes normal 


Subroutine listing: 


SUBROUTINE PDIV <P, IDIMF,X,IDIMX,Y,IDIMY,TOL* IER) 
C THIS ROUTINE OIVIDES TWO POLYNOMIALS RETURNING 

C THE QUOTIENT ANO THE REMAINDER 

0 IMENS ION P(l) ,X( 1) ,Y( 1) 

5 10 I DI HP* ID I MX-I0IMY+1 

IF( IDIMP >20*30 *60 
20 1 01 MP= 0 
30 IER=0 
40 RETURN 

10 00 IER=1 

GO TO_40j 
60 ioi'm'x=ioimy-i 

1= I DI M F 
70 II* I+I DIMX 

15 P(I)*X (III /YtlDIMY) 

00 80 K=1 » IOIMX 
J=K-1+I 

80 X(J>=X<J)-P<I>*YOO 
1 = 1-1 

20 IFC 1)90,90*70 

90 CALL PNORMCX,IDIMX*TCL) 

GO TO 30 
END 


Subroutine PNORM (X, IDIMX, EPS) 


Purpose: To normalize coefficients of a polynomial 


Subroutine arguments: 


X array of coefficients for p (oc) 

IDIMX dimension of X; replaced by final dimension after normalization 

EPS tolerance below which coefficient is eliminated 


Subroutine listing: 


SUBROUTINE PNORMt X* IDIMX, EPS) 
DIMENSION X(l> 

1 IF<I0IMX>4,4,2 

2 IF< ABS(X(I0IMX))-EPS) 3*3,4 

5 3 IDIMX* ID IMX-1 

GO TO 1 
4 RETURN 
END 
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APPENDIX - Continued 


Subroutine ROOE (A , B , N , BOUND , ROOT , NI) 
Purpose: To compute all the real zeros of p Or) 

Flow chart: 



Subroutine arguments: 

A array of coefficients for p ( x ) 

B work array 

N degree of p (x) 

BOUND closed interval [-BOUND .BOUND] which contains all the real 
zeros of p (x) , an arbitrary polynomial 
ROOT array containing the real zeros 

NI number of real zeros in array root 
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APPENDIX - Continued 


Subroutine listing: 


5 


10 


15 


20 


25 


30 


35 


40 


45 


50 


55 


C 

C 

c 

c 

c 

c 

c 


/ 

SUBROUTINE R00E ( A • B ,N. BOUND , RC CT, ft I 1 
DIMENSION e<i),A(l >,ROOT(l) 

DIMENSION Y (30 ) 

REAL LX 



THIS routine actually COMPUTES THE real ROOTS. 

e,dr,xl, are tolerances upon which convergence is based;their 

ORDER INDICATING THEIR SEVERITY, MORE THAN ONE IS USED TO IMPROVE 
PROBABILITY OF CONVERGENCE WHEN ACCURACY IS LOST AND/OR PRESENCE 
OF MULTIPLE ROOTS. 


JJ=N*10 
KB=NI * 
/ E=. 1**16 
' GK=0 « 



0 = 0 . If 
XL=.1**6 • 
P=0 

G=. 1**10 
OR=Q • 

N 1= 0 
LK-N+1 » 
LL=LK 


DO 190 J=1 , LK 
190 B ( J 1=A (J l 
XX= 1/B (11 


vv=xx 

DO 444 J=1 * LK 
444 B (J 1 =B ( J 1 * XX 
1 = 0 

183 L J=2 

I J=K0-I 

ANS=0 

G=0 

IF( I J. NE. 11G0 TO 2566 
JH^ALOGIO (EOUNC)/ALCG(2. 1 
IN=JH/6 
C 

C IF ONLY ONE DISTINCT ZERO REMAINS, THEN ISCLATE IT. 

C 

DO 117 J=l,10j)Q- 
G=B0UND-(BCUNQ/(2.**(J*IM) ) 1 
SS=BOUND-G- 

IF(SS.LT.100. > GO TO 1000 
CALL CHECK (B,LK-1,G,IER) 

I F ( I ER *NE • 0 ) GO TO 1000 
117 CONTINUE 
1000 CONTINUE 

G=BOUN0- (BOUND/ (2.** C ( J-1)*IM) ) ) 

2566 CONTINUE 
XrSOUNO-G 
J 182 CONTINUE 
C 

C DETERMINES STEP SIZE ,DEG. 

C 

0EG=(2*X>/ (2.**(LJ-11 1 
LG=2*X/DEG 
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APPENDIX - Continued 


60 


65 


70 


75 


80 


85 


90 


95 


100 


105 


110 


T = 0 
150 1=2 

22 LX=M00(L*2>-1. 

C 

C COMPUTES INITIAL GUESSES WITH STEP SIZEtDEG. 

C 

XO*SIGMG.LX)-(SIGMDEG*LX)»( (L/2)-i> ) 
Xl*SIGN(GtLX)-(SIGN(0EG*LX)ML/2l I 
C 

CALL PVAL (Z*X1,B,LK) 

Z=ABS( Z) +10. 

CALL PVAL(R.X0,8,LK> 

20 CONTINUE 
C 

C EVALUATES F<X> AT COMPUTEO ITERATE 

C 

CALL PVAL< S* X1»B*LKI 
Iff IJ.NE.IJGO TO 103 
II=K8-I 

IF ( II. GT • 1 IGO TO 103 
F=SIGN(1. ,S> 

H=SIGN(1. ,RI 
IF (F*H.GT. CIGO TO 103 
XM- (S-R I / ( Xl-XO) 

CT=A0S(X1-XO> 

IF ( CT • GT. . 5)GO TO 97 
IF(CT.GT..01*ABSCX1))GO TO 97 
XMM=S-XM*X1 
ANS=-XMM/XM 
GC TO 206 
97 CONTINUE 
103 CONTINUE 
C 

C DETERMINES MINIMUM F(X) IF DIVERGENCE OCCURS FOR ALL ITERATES. 

C 

IF(T.EQ.Q)GO TO 122 
I F ( A8S (SI.GT.TJGC TO 10 
122 CONTINUE 
T=A8S(S) 

G 1= X 1 

10 CONTINUE 
C 

C IF F ( X > CONVERGES . THEN ISOLATE ZERO ANO REDUCE P<XI 

C 

IFCABS(S) • LT. E ) GO TO 125 

IF((A8S(S) .GT.Z>. AND. (A9S(S) .LT.DR) IGO TO 125 
IFCABS (SI .LT,XL)f>=Xl 
C 

C CHECKING FOR DIVERGENCE 

C 

I F ( ABS (S I • GT. 2 1 GO TO 120 

Z=ABS (SI 

V=Xi 

IF(S-R.EQ.O)GO TO 120 
FG=(X1-X0) /(S-R) 

X 1=X1- (F G*S) 

IF (ABS (XI) .GT.80UN0)X1=SIGN(BCUNC,X1> 



APPENDIX — Continued 


■III llllll III I 


115 


120 


125 


130 


135 


140 


145 


150 


155 


160 


165 


170 


175 


160 


185 


190 


195 


IF CABS (XI) .LT.G)X1=SIGMG,X1> 

xo=v 

GO S TO 20 
125 1=1+1 
ANS=Q 
P=Q 
NI=I 

23 CONTINUE 
ROOT C I )=X1 
IFCI.EQ.NIGO TC 206 
IFCGK.NE.O.JGO TO 105 
C 

C COMPUTES OEFLATEO POLY NOMI AL ,USING SYNTHETIC OIVISION, DETERMINE 

C NEW INTERVAL. 

C 

550 CONTINUE 

CALL S YNTH <8* LK ♦ XI ) 

IF (LK-l*EQ#2)GO TO 100 

LK=LK-1 

C 

C OETERMINE 80UNG FOR ZERO ANO HENCE THE NO. OF OISTINCT ZEROES 

C FOR THE RECUCEO FOLYHOMIAL 

C 

CALL BOUNE C8, LK-1 » BOUND , Y) 

CALL CHECK(B,LK-1,80UN0,IER) 

C 

C interrogates THE STATUS OF COMPLETION, 

c 

IFCI.LT. KB)GO TO 2300 

IF< (IER.EQ.O) .AND. (I.GE.KB) >GO TO 206 
2300 CONTINUE 
V V= i/B Cl) 

DO 355 J=1,LK 
355 8CJ>=8CJ)*VV 
699 CONTINUE 
GO TO 18 3 
100 Xl=-BC 1) /B C2 ) 

GO TO 125 
120 L=L+1 
C 

C COMPUTES THE NEXT INITIAL GUESS IF INTERVAL IS NOT EXHAUSTEO. 

C 

IF(L.LT.LG)G0 TO 22 
• IX=IX+1 

IF CP* EO. 0 ) GO TO 1212 
X1=P 

GO TO 125 
1212 CONTINUE 
C 

C IF STEP SI2E IS SMALL ENOUGH, SUBJECT POLYNOMIAL TO BISECTION, 

C OTHERWISE, CHOOSE A SMALLER STEP SIZE, ETC, 

C 

IFCDEG.LT. ,1>G0 TO 1556 
LJ=LJ+i 
GO TO 182 
1556 CONTINUE 
C 

C 8EGIN BISECTION ARGUMENT BASED ON Gi AND DEG. 

C 

X i=Gl 
D=D£G' 

JH=ALQG10CCEG) /ALOG10C2.) 

JH= JH+18 
DO 18 J=1,JH 
XC=Xl+0 
XG=Xl-0 

CALL PVALCS,XC,A,LLI 
S=ABS(S1 

CALL PVAL CT ,XG, A, LL) 

T=A 8S ( Ti 

Xl=Xl-CSIGMD/2» S-TI ) 

CALL PVAL C H, XI , A, LL) 

H*A8S C H) 

IFCH.LT.XUGO TO 125 
0 = 0/2 

18 CONTINUE 
206 CONTINUE 

IF ( ANS«EO. 0)GO TO 105 

Xi=ANS 

GK=1 

GO TO 125 
105 CONTINUE 
RETURN 
END 
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APPENDIX - Concluded 


Subroutine SYNTH (B , K , X) 


Purpose: To deflate a polynomial by one degree using synthetic division 


Subroutine arguments: 


B array of coefficients for a polynomial , p (x) 

K degree + 1 

X zero of p (oe) 

Subroutine listing: 


s 


10 


SUBROUTINE SYNTH<8,K,X1 
(DIMENSION B <1 ) «C ( 30 ) 

C 

c this algorithm oeelates the polynomial using synthetic division. 
c 

L=K-1 

C(K)=Q(K> 

0C 17 J=1*L 

17 C <K-J1=X*C (K-J«-l> *8<K-J> 

00 18 J=1,L 

18 0(J)=C(J-H) 

return 

ENO 
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